﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Bouyei.GeoCore.Geometries
{
    using GeoParsers;

    public class GeoEllipsoid : GeoBase
    {
        readonly double ep = 0;
        readonly double efourp = 0;//pow(4)
        readonly double esixp = 0;//pow(6)
        readonly double eeightp = 0;//pow(8)

        readonly double A = 0, B = 0, C = 0, D = 0, E = 0;

        public GeoEllipsoid(int zoneWide = 3, int precision = 0)
       : base(zoneWide: zoneWide, precision: precision)
        {
            ep = ((lRadius * lRadius - sRadius * sRadius) / (lRadius * lRadius));
            efourp = ep * ep;// ep * ep;
            esixp = efourp * ep;// efourp * ep;
            eeightp = efourp * efourp;// efourp * efourp;
            //sRadiusp = (sRadius * sRadius);

            //A = 1 + 0.5 * ep + (30 / 80.0) * efourp + (35 / 112.0) * esixp + (630 / 2304.0) * eeightp;
            //B = (1 / 6.0) * ep + (15 / 80.0) * efourp + (21 / 112.0) * esixp + (420 / 2304.0) * eeightp;
            //C = (3 / 80.0) * efourp + (7 / 112.0) * esixp + (180 / 2304.0) * eeightp;
            //D = (1 / 112.0) * esixp + (45 / 2304.0) * eeightp;
            //E = (5 / 2304.0) * eeightp;

            //简化
            A = 1 + 0.5 * ep + 0.375 * efourp + 0.3125 * esixp + 0.2734375 * eeightp;
            B = 0.16666666666667 * ep + 0.1875 * efourp + 0.1875 * esixp + 0.18229166666667 * eeightp;
            C = 0.0375 * efourp + 0.0625 * esixp + 0.078125 * eeightp;
            D = 0.00892857142857 * esixp + 0.01953125 * eeightp;
            E = 0.00217013888889 * eeightp;
        }

        public GeoEllipsoid(double lRadius, double sRadius, int zoneWide = 3, int precision = 0)
            : base(lRadius, sRadius, precision, zoneWide)
        {
            ep = ((lRadius * lRadius - sRadius * sRadius) / (lRadius * lRadius));
            efourp = ep * ep;// ep * ep;
            esixp = efourp * ep;// efourp * ep;
            eeightp = efourp * efourp;// efourp * efourp;
            //sRadiusp = (sRadius * sRadius);

            //A = 1 + 0.5 * ep + (30 / 80.0) * efourp + (35 / 112.0) * esixp + (630 / 2304.0) * eeightp;
            //B = (1 / 6.0) * ep + (15 / 80.0) * efourp + (21 / 112.0) * esixp + (420 / 2304.0) * eeightp;
            //C = (3 / 80.0) * efourp + (7 / 112.0) * esixp + (180 / 2304.0) * eeightp;
            //D = (1 / 112.0) * esixp + (45 / 2304.0) * eeightp;
            //E = (5 / 2304.0) * eeightp;

            //简化
            A = 1 + 0.5 * ep + 0.375 * efourp + 0.3125 * esixp + 0.2734375 * eeightp;
            B = 0.16666666666667 * ep + 0.1875 * efourp + 0.1875 * esixp + 0.18229166666667 * eeightp;
            C = 0.0375 * efourp + 0.0625 * esixp + 0.078125 * eeightp;
            D = 0.00892857142857 * esixp + 0.01953125 * eeightp;
            E = 0.00217013888889 * eeightp;
        }

        public double Area(Coordinate[] coords)
        {
            double s = 0;
            int len = coords.Length - 1;

            bool isXY = coords[0].X > 360;

            if (isXY)
            {
                for (int i = 0; i < coords.Length; ++i)
                {
                    coords[i] = XYtoLB(coords[i]);
                }
            }

            for (int i = 0; i < len; ++i)
            {
                s += Trapezoid(coords[i], coords[i + 1]);
            }

            if (precision > 0)
            {
                s = Math.Round(s, precision);
            }

            return s;
        }

        public double Area(string wkt)
        {
            GeoParsers.Wkt.GeometryReader wktReader = new GeoParsers.Wkt.GeometryReader(wkt);
            var colleciton = wktReader.Reader();

            var sum =Area(colleciton[0]);
            double holes = 0;
            //图形有洞
            for (int i = 1; i < colleciton.Count; ++i)
            {
                holes +=Area(colleciton[i]);
            }
            return Math.Abs(sum + holes);
        }

        public double Area(Geometry geometry)
        {
            double sum = 0;

            for (int i = 0; i < geometry.GeometryCount; ++i)
            {
                var geo = geometry.GetGemoetry(i);
                for (int j = 0; j < geo.Count; ++j)
                {
                    var seq = geo[i];
                    sum += Area(seq.GetCoordinates());
                }
            }

            return Math.Abs(sum);
        }

        /// <summary>
        /// 计算梯形面积
        /// </summary>
        /// <param name="start"></param>
        /// <param name="end"></param>
        /// <returns></returns>
        private double Trapezoid(Coordinate start, Coordinate end)
        {
            double x = ToRadian(start.X);
            double y = ToRadian(start.Y);
            double x1 = ToRadian(end.X);
            double y1 = ToRadian(end.Y);

            double L = (x1 + x) * 0.5;//文档说明是经度差? 这里使用经度平均值
            double Lm = (y1 - y) * 0.5;//纬差
            double Bm = (y1 + y) * 0.5;//维度平均值

           double  s = 2 * sRadius * L*sRadius * (
                A * Math.Sin(Lm) * Math.Cos(Bm)
                - B * Math.Sin(3 * Lm) * Math.Cos(3 * Bm)
                + C * Math.Sin(5 * Lm) * Math.Cos(5 * Bm)
                - D * Math.Sin(7 * Lm) * Math.Cos(7 * Bm)
                + E * Math.Sin(9 * Lm) * Math.Cos(9 * Bm));

            return s;
        }
    }
}
